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Abstract 

We perform lattice simulations of two-flavor QCD using Neuberger's overlap fermion, with which 
the exact chiral symmetry is realized at finite lattice spacings. The e-regime is reached by decreasing 
the light quark mass down to 3 MeV on a 16'^ x 32 lattice with a lattice spacing ~ 0.11 fm. We 
find a good agreement of the low-lying Dirac eigenvalue spectrum with the analytical predictions of 
the chiral random matrix theory, which reduces to the chiral perturbation theory in the e-regime. 
The chiral condensate is extracted as S^^(2 GeV) = (251 ± 7 ± 11 MeV)^, where the errors are 
statistical and an estimate of the higher order effects in the e-expansion. 
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In Quantum Chromodynamics (QCD), it is widely believed that chiral symmetry is spon- 
taneously broken, making pions nearly massless while giving masses of order Aqcd, the QCD 
scale, to the other hadrons. In fact, the Chiral Perturbation Theory (ChPT), an effective 
theory based on the spontaneously broken chiral symmetry, describes low energy interac- 
tions of pions very accurately. Nevertheless, theorists have not been successful in analytically 
solving QCD and deriving the chiral symmetry breaking, due to its highly non-perturbative 
dynamics. 

The most promising approach to establishing the link between QCD and ChPT is to 
utilize the numerical simulation of lattice QCD, with which the every prediction of ChPT 
can be tested in principle. For instance, the presence of the so-called chiral logarithms, the 
effect of pion cloud, should be reproduced. Such a numerical test is, however, not an easy 
task, because of rapidly increasing computational cost in the small quark mass region where 
ChPT is reliably applied. Another serious problem is the explicit violation of the chiral 
symmetry at finite lattice spacings in the conventional fermion formulations, with which 
the conclusive test of ChPT requires well-controlled and thus computationally demanding 
continuum extrapolation. 

In this work we improve this situation in two ways. First, we employ Neuberger's overlap 
fermion for dynamical quarks. It preserves exact chiral symmetry at finite lattice 

spacings, and hence ChPT can be applied before taking the continuum limit. Although 
the numerical cost of the overlap fermion is almost 100 times higher than that of the other 
fermions, new computational facilities at KEK enable us to carry out such a work. 

Second, we study the correspondence between QCD and ChPT in the so-called e-regime 

0.QQ 

, which is characterized by the small pion mass satisfying rriT^L ^ 1 with L the 
box size. In this regime ChPT is safely applied as an expansion in terms of ~ itIt^/Aqcd, 
provided that the condition 1/(Aqcd-^)^ ^ 1, the usual condition that the box size is larger 
than the inverse QCD scale, is satisfied. We set the sea quark mass to ~ 3 MeV, for which 
m^L ~ 1.0. With the space-time volume x T ^ (1.8 fm)^ x (3.5 fm), the numerical cost is 
still not prohibitive even with such a small sea quark mass, since the finite volume provides 
a natural lower bound on the lowest eigenvalue of the Dirac operator. 

In the e-regime, zero-momentum modes of the pion field dominate the dynamics and 
the kinetic term gives only subleading contributions. The Lagrangian of ChPT reduces to 
£(0) _ j77,xi ReTr[f/] with m the (degenerate) quark mass and U (g SU(A^/)) the pion field 
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{Nf = 2 in this work). The system is fully characterized by a parameter mHV, where S is 
the chiral condensate and V is the space-time volume x T. Dependence on the topological 
charge Q of the gauge field also becomes significant. 

At the leading order of the e-expansion, the chiral Random Matrix Theory (ChRMT) 
provides an equivalent description of ChPT Furthermore, ChRMT can predict 

the distributions of the individual eigenvalues of the Dirac operator, which may be directly 
compared with the lattice data. In the quenched approximation, a good agreement between 
ChRMT and lattice calculation has been observed using the overlap-Dirac operator 



HI], and S has been determined by matching the eigenvalues 12[. The present work is an 



extension of these works to two-flavor QCD. A preliminary report of this work has been 
presented in [l^, and an overview of our dynamical overlap fermion project is found in fl^ . 



but the e-regime was not reached because 



Similar studies have been done recently 
of larger sea quark masses. 

We have performed numerical simulations on a 16^ x 32 lattice at a lattice spacing a 



0.11 fm as determined 



the overlap fermion 



rom the scale ro(=0.49fm) of the heavy quark potential. We employ 
21], whose Dirac operator is 



D{m) 



m 

mo + - 



rriQ 



Tfl \ 

— \ 75Sgn[if,4/(-mo)] 



(1) 



for the quark mass m. Here, Hw{—niQ) denotes the standard hermitian Wilson-Dirac op- 
erator Hy^{—mQ) = '-^^D^Y^—mo) with a large negative mass term (we choose uiq = 1.6 
throughout this work). For the gauge part, the Iwasaki action is used at P = 2.35 together 



171 ]. which preserves 



with unphysical Wilson fermions and associated twisted-mass ghosts 
the global topological charge during molecular-dynamics evolutions of the gauge field. This 
is desirable for the e-regime simulations since we can effectively accumulate statistics at a 
given topological charge. In this work our simulation is confined in a topological sector 
Q = 0. 



For the simulation with the dynamical overlap fermions |18|, we use the Hybrid Monte 
Carlo (HMC) algorithm. The sign function in ([1]) is approximated by a rational function with 
Zolotarev's optimal coefficients after projecting out low-lying eigenvalues of ]ifvK(— mo)]. 
With 10 poles the sign function has a lO^^"^^*) precision. Thanks to the extra Wilson 
fermions, the lowest eigenvalue of Hwi—iTio) never passes zero, and hence no special care of 
the discontinuity of the fermion determinant is needed. 



3 



b 6x10 



I 4x10 




FIG. 1: Number of the Wilson-Dirac operator multiplication per trajectory (upper panel) and 
per an overlap inversion (lower panel). The curves are fit to data above m = 0.030 with the form 
oc l/m". 



The simulation cost is substantially reduced by the mass preconditioning of the HMC 
Hamiltonian j2o|. The heavier overlap fermion mass for the preconditioner is chosen to 0.4 
except for two lightest sea quark masses where the value is 0.2. The relaxed conjugate 
gradient algorithm to invert the overlap-Dirac operator 19|] also helps to speed up the 
simulation by about a factor of 2. 

For the sea quark mass we take seven values: 0.110, 0.090, 0.065, 0.045, 0.030, 0.020, and 
0.002. The lightest sea quark (m = 0.002) corresponds to the e-regime. The simulation cost 
measured by the number of the Wilson-Dirac operator multiplication is plotted in Fig. [H 
The upper panel shows the cost per trajectory of length 0.5; the lower panel presents the 
cost of inverting the overlap-Dirac operator when we calculate the Hamiltonian at the end 
of each trajectory. Increase of the numerical cost towards the chiral limit is not as strong 
as expected: if we fit the data assuming the scaling law ~ 1/m" above m = 0.030, the 
power a is 0.82 for the inversion and 0.49 for the total cost. In the e-regime the cost is even 
lower than the expectation from the power law. This is because the cost is governed by the 
lowest-lying eigenvalue rather than the quark mass as explained below. 

The number of trajectories is 1,100 for each sea quark mass after discarding 500 for 
thermalization. For the e-regime run at m = 0.002, we have accumulated 4,600 trajectories. 
In addition, we have generated quenched lattices at a similar lattice spacing with Q = and 
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FIG. 2: Lowest 5 eigenvalues as a function of sea quark mass. Dashed line shows X = m. 

2. The computational cost at m = 0.002 is about one hour per trajectory on a half rack 
(512 nodes) of the IBM BlueGene/L (2.8 TFlops peak performance). 

The lowest 50 eigenvalues of the overlap-Dirac operator -D(O) are calculated at every 
10 trajectories. We use the implicitly restarted Lanczos algorithm for a chirally projected 
operator -D(O) P+, where P_|_ = (l+75)/2. From its eigenvalue ReA"*" the pair of eigenvalues 
X""" (and its complex conjugate) of -D(O) is extracted through the relation |1 — A°^/mop = 1, 
that forms a circle on a complex plane. For the comparison with CliRMT, the lattice 
eigenvalue X°'" is projected onto the imaginary axis as A = ImA°''/(l — ReA°^/(2mo)). Note 
that A is very close to ImA°^ (within 0.05%) for the low-lying modes we are interested in. 

In Fig. [2] we plot the ensemble averages of the lowest 5 eigenvalues (A^) {k = 1-5) as 
a function of the sea quark mass. We observe that the low-lying spectrum is lifted as 
the chiral limit is approached. This is a direct consequence of the fermion determinant 
~ rifcdAfcP + m^), which repels the small eigenvalues from zero when m becomes smaller 
than the lowest eigenvalue. This is exactly the region where the numerical cost saturates as it 
is controlled by Ai rather than m. We also find that the autocorrelation length of the lowest 
eigenvalues is significantly longer in the e-regime. We therefore use the jackknife method 
in the statistical analysis with a bin size of 200 trajectories, with which the statistical error 
saturates. 

In the e-regime, ChRMT predicts the probability distribution Pk{Ck) of lowest-lying eigen- 
values (k, and thus their ensemble averages {(k) = /o°° '^CfcCfcPfc(Cfc)- The correspondence 
between ChRMT and ChPT implies the relation {(k) = (Afc)SV", from which we can extract 
S, once (Afc)'s are obtained. Correction for finite sea quark mass m can also be calculated 
by ChRMT. At m = 0.002, it is only 0.9% for the lowest eigenvalue, which is taken into 
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FIG. 3: Ratio of the eigenvalues {Ck) / iCi) ^or combinations of k and I G 1-4 (denoted in the plot 
as k/l). In addition to the two-flavor QCD data (middle), quenched data at Q = (left) and 2 
(right) are also shown. Lattice data (circles) are compared with the ChRMT predictions (bars). 

account in the following analysis. 

We first compare the pattern of the eigenvalue spectrum of the Dirac operator. In the 
ratios {(k) / (0) of low-lying eigenvalues the factor J^V drops off and the comparison is pa- 
rameter free. As Fig. [3] shows, the lattice data agree well with the ChRMT predictions 
(middle panel). It is also known that there exists the so-called flavor-topology duality in 
ChRMT: the low-mode spectrum is identical between the two-flavor (massless) theory at 
Q = and the quenched theory at Q = 2 (right panel), while the quenched spectrum at 
Q = is drastically different (left panel). This is nicely reproduced by the lattice data. 
These clear patterns indicate that the low-lying modes of the QCD Dirac operator are in 
fact responsible for the zero-momentum mode of the pion field U, which induce the vacuum 
with spontaneously broken chiral symmetry. 

If we look into the details, however, there is a slight disagreement in higher eigenvalues. 
For instance, the deviation in the combination k/l = 4/1 is 10%. The reason for this will 
be discussed later. 

Another non-trivial comparison can be made through the shape of the eigenvalue dis- 
tributions. We plot the cumulative distribution Cfc(Cfc) = Jo'' ^CpkiC) of the three lowest 
eigenvalues in Fig. HI The agreement between the lattice data and ChRMT (solid curves) 
is quite good for the lowest eigenvalue. For the higher modes the agreement of the central 
value is marginal while the shape seems well described by ChRMT. Table H] lists the numer- 
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FIG. 4: Cumulative distribution of the low- lying eigenvalues. The horizontal error comes from 
the statistical error of S. The solid curves are the ChRMT results with the input from the average 
(Ai). 









k 


= 1 4.30 [4.30] 


1.234 1.215(48) 


k 


= 2 7.62 7.25(13) 


1.316 1.453(83) 


k 


= 3 10.83 9.88(21) 


1.373 1.587(97) 


k 


= 4 14.01 12.58(28) 


1.414 1.54(10) 



TABLE I: Comparison of the low-mode spectrum of ChRMT and the lattice data. The number 
given in [] is used as an input. 

ical results of both ChRMT and lattice data. To characterize the shape of the distribution, 
we define the width (A^^) = \/ id) " (Cfe)^ well as its lattice counterpart {A\k)T.V. The 
overall agreement is very good, but we see deviations of about 10% in the averages and 16% 
in the widths, at the largest. 

Because of the exact chiral symmetry in our simulation, we do not expect significant 
effects due to finite a in the comparison of low-lying eigenvalues. The largest possible source 
of systematic errors is the higher order effects in the e-expansion, which is a finite volume 
effect. Although the higher order corrections can not be calculated within the framework 
of ChRMT, an estimate of their size can be given by ChPT. The next-to-leading order 
correction to the chiral condensate is given by S[l -|- {{Nf — 1)/Nf)[3i/{FLY], where [3i is 
a numerical constant depending on the lattice geometry [Sj. Numerically, the correction is 
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13% (we assume the pion decay constant F = 93 MeV), which is about the same size of the 
deviation of the eigenvalue distributions. 

Since the contribution from finite momentum modes is expected to be more significant 
for higher eigenvalues, we take the lowest eigenvalue as an input for the determination of 
S. From the average of Ai we obtain S'"*a^ = 0.00212(6) in the lattice unit and S''** = 
[240(2) (6) MeVf in the physical unit. The second error in the latter comes from the lattice 
scale a = 0.107(3) fm. We put a superscript lat in order to emphasize that it is defined on 
the lattice. 

In order to convert S'"* to the definition in the MS scheme, we have calculated the 
renormalization factor Z^^{2 GeV) using the non-perturbative renormalization technique 
through the RI/MOM scheme [21|]. Calculation is done at m = 0.002 with several different 
valence quark masses. The result is Zf^{2 GeV) = 1.14(2). Details of this calculation will 
be presented elsewhere. 

Including the renormalization factor, our final result is S'^^(2 GeV) = [251(7) (11) MeV]^. 
The errors represent a combined statistical error (from Ai, tq and Zg^^) and the systematic 
error estimated from the higher order effects in the e-expansion, respectively. Since the 
calculation is done at a single lattice spacing, the discretization error cannot be quantified 
reliably, but we do not expect much larger error because our lattice action is free from 0(a) 
discretization effects. 

In this letter we demonstrated that the lattice QCD simulation is feasible near the chiral 
limit, as far as the exact chiral symmetry is respected. The link between QCD and ChPT is 
established in the e-regime without recourse to chiral extrapolations. From their correspon- 
dence, the low energy constants, such as the chiral condensate, can be precisely calculated. 
Improvement of the precision may be achieved by increasing the physical volume, as the 
higher order effects in the e-expansion are suppressed as 1/L^. Further information on the 
low energy constants can be extracted in the e-regime by calculating two- and three-point 
unctions or analyzing the Dirac eigenvalue spectrum with imaginary chemical potential 



22, l23|, l2J| . This work is a first step towards such programs. 
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